home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_Jnu.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.0 KB  |  178 lines
  1. /* specfunc/bessel_Jnu.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_bessel.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "bessel.h"
  30. #include "bessel_olver.h"
  31. #include "bessel_temme.h"
  32.  
  33.  
  34. /* Evaluate at large enough nu to apply asymptotic
  35.  * results and apply backward recurrence.
  36.  */
  37. #if 0
  38. static
  39. int
  40. bessel_J_recur_asymp(const double nu, const double x,
  41.                      gsl_sf_result * Jnu, gsl_sf_result * Jnup1)
  42. {
  43.   const double nu_cut = 25.0;
  44.   int n;
  45.   int steps = ceil(nu_cut - nu) + 1;
  46.  
  47.   gsl_sf_result r_Jnp1;
  48.   gsl_sf_result r_Jn;
  49.   int stat_O1 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps + 1.0, x, &r_Jnp1);
  50.   int stat_O2 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps,       x, &r_Jn);
  51.   double r_fe = fabs(r_Jnp1.err/r_Jnp1.val) + fabs(r_Jn.err/r_Jn.val);
  52.   double Jnp1 = r_Jnp1.val;
  53.   double Jn   = r_Jn.val;
  54.   double Jnm1;
  55.   double Jnp1_save;
  56.  
  57.   for(n=steps; n>0; n--) {
  58.     Jnm1 = 2.0*(nu+n)/x * Jn - Jnp1;
  59.     Jnp1 = Jn;
  60.     Jnp1_save = Jn;
  61.     Jn   = Jnm1;
  62.   }
  63.  
  64.   Jnu->val = Jn;
  65.   Jnu->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jn);
  66.   Jnup1->val = Jnp1_save;
  67.   Jnup1->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jnp1_save);
  68.  
  69.   return GSL_ERROR_SELECT_2(stat_O1, stat_O2);
  70. }
  71. #endif
  72.  
  73.  
  74. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  75.  
  76. int
  77. gsl_sf_bessel_Jnu_e(const double nu, const double x, gsl_sf_result * result)
  78. {
  79.   /* CHECK_POINTER(result) */
  80.  
  81.   if(x < 0.0 || nu < 0.0) {
  82.     DOMAIN_ERROR(result);
  83.   }
  84.   else if(x == 0.0) {
  85.     if(nu == 0.0) {
  86.       result->val = 1.0;
  87.       result->err = 0.0;
  88.     }
  89.     else {
  90.       result->val = 0.0;
  91.       result->err = 0.0;
  92.     }
  93.     return GSL_SUCCESS;
  94.   }
  95.   else if(x*x < 10.0*(nu+1.0)) {
  96.     return gsl_sf_bessel_IJ_taylor_e(nu, x, -1, 100, GSL_DBL_EPSILON, result);
  97.   }
  98.   else if(nu > 50.0) {
  99.     return gsl_sf_bessel_Jnu_asymp_Olver_e(nu, x, result);
  100.   }
  101.   else {
  102.     /* -1/2 <= mu <= 1/2 */
  103.     int N = (int)(nu + 0.5);
  104.     double mu = nu - N;
  105.  
  106.     /* Determine the J ratio at nu.
  107.      */
  108.     double Jnup1_Jnu;
  109.     double sgn_Jnu;
  110.     const int stat_CF1 = gsl_sf_bessel_J_CF1(nu, x, &Jnup1_Jnu, &sgn_Jnu);
  111.  
  112.     if(x < 2.0) {
  113.       /* Determine Y_mu, Y_mup1 directly and recurse forward to nu.
  114.        * Then use the CF1 information to solve for J_nu and J_nup1.
  115.        */
  116.       gsl_sf_result Y_mu, Y_mup1;
  117.       const int stat_mu = gsl_sf_bessel_Y_temme(mu, x, &Y_mu, &Y_mup1);
  118.       
  119.       double Ynm1 = Y_mu.val;
  120.       double Yn   = Y_mup1.val;
  121.       double Ynp1 = 0.0;
  122.       int n;
  123.       for(n=1; n<N; n++) {
  124.         Ynp1 = 2.0*(mu+n)/x * Yn - Ynm1;
  125.     Ynm1 = Yn;
  126.     Yn   = Ynp1;
  127.       }
  128.  
  129.       result->val = 2.0/(M_PI*x) / (Jnup1_Jnu*Yn - Ynp1);
  130.       result->err = GSL_DBL_EPSILON * (N + 2.0) * fabs(result->val);
  131.       return GSL_ERROR_SELECT_2(stat_mu, stat_CF1);
  132.     }
  133.     else {
  134.       /* Recurse backward from nu to mu, determining the J ratio
  135.        * at mu. Use this together with a Steed method CF2 to
  136.        * determine the actual J_mu, and thus obtain the normalization.
  137.        */
  138.       double Jmu;
  139.       double Jmup1_Jmu;
  140.       double sgn_Jmu;
  141.       double Jmuprime_Jmu;
  142.       double P, Q;
  143.       const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);
  144.       double gamma;
  145.  
  146.       double Jnp1 = sgn_Jnu * GSL_SQRT_DBL_MIN * Jnup1_Jnu;
  147.       double Jn   = sgn_Jnu * GSL_SQRT_DBL_MIN;
  148.       double Jnm1;
  149.       int n;
  150.       for(n=N; n>0; n--) {
  151.         Jnm1 = 2.0*(mu+n)/x * Jn - Jnp1;
  152.         Jnp1 = Jn;
  153.         Jn   = Jnm1;
  154.       }
  155.       Jmup1_Jmu = Jnp1/Jn;
  156.       sgn_Jmu   = GSL_SIGN(Jn);
  157.       Jmuprime_Jmu = mu/x - Jmup1_Jmu;
  158.  
  159.       gamma = (P - Jmuprime_Jmu)/Q;
  160.       Jmu   = sgn_Jmu * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jmuprime_Jmu)));
  161.  
  162.       result->val = Jmu * (sgn_Jnu * GSL_SQRT_DBL_MIN) / Jn;
  163.       result->err = 2.0 * GSL_DBL_EPSILON * (N + 2.0) * fabs(result->val);
  164.  
  165.       return GSL_ERROR_SELECT_2(stat_CF2, stat_CF1);
  166.     }
  167.   }
  168. }
  169.  
  170. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  171.  
  172. #include "eval.h"
  173.  
  174. double gsl_sf_bessel_Jnu(const double nu, const double x)
  175. {
  176.   EVAL_RESULT(gsl_sf_bessel_Jnu_e(nu, x, &result));
  177. }
  178.